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Abstract 

Many biological, physical, and social interactions have a particular dependence on where they take 
place. In living cells, protein movement between the nucleus and cytoplasm affects cellular response 
(i.e., proteins must be present in the nucleus to regulate their target genes). Here we use recent devel- 
opments from dynamical systems and chemical reaction network theory to identify and characterize 
the key-role of the spatial organization of eukaryotic cells in cellular information processing. In par- 
ticular the existence of distinct compartments plays a pivotal role in whether a system is capable of 
multistationarity (multiple response states), and is thus directly linked to the amount of information 
that the signaling molecules can represent in the nucleus. Multistationarity provides a mechanism for 
switching between different response states in cell signaling systems and enables multiple outcomes 
for cellular-decision making. We find that introducing species localization can alter the capacity for 
multistationarity and mathematically demonstrate that shuttling confers flexibility for and greater con- 
trol of the emergence of an all-or-none response. 

Keywords: chemical reaction networks, MAPK, mass-action kinetics, cellular information pro- 
cessing, spatial localization 

Introduction 

Cells constantly have to adapt and respond to their environment. In single-celled organisms those cells 
least well adjusted to their surroundings will tend to contribute less to future generations than cells that are 
able to assimilate better or more quickly to changing circumstances. In multi-cellular organisms, aberrant 
response of individual cells to environmental or physiological cues may result in developmental anomalies 
or disease. The way in which cells respond to external signals, process them and act upon them is thus 
intimately and inextricably linked to an organism's fate (HE); and in the long-term evolution will shape 
the molecular machinery underlying cellular decision making processes J3HU. 

One central aspect of biological information processing is the mapping of environments onto intra- 
cellular states given by the abundances of the molecular species (proteins, mRNAs, metabolites etc) under 
consideration. In this processing of information one or more environmental variables need to be repre- 
sented in a way that facilitates the appropriate response. Continuous and discrete representations have 
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been reported, and it is easy to see how an increased number of states will start to mimic the "analog" 
nature of continuously varying states; here, of course, we are typically only interested in stable states. A 
simple "on/off switch, for example, is a binary representation or response mechanism; this behavior is 
particularly interesting if there is a regime of conditions where the system can populate either state. In this 
case we speak of a bistable switch; outside this regime we only find a single state for the system. More 
generally we speak of multistability if more than two stable states are obtainable simultaneously. 

The number of states in which a cell can be at any given time is linked to the flexibility in its decision 
making [5 |. If only one state is accessible (and stable) then there is obviously no room for "choice" (even 
if such choices are made by random processes) and any cell-to-cell variability will derive from intrinsic or 
extrinsic sources of noise which will broaden out the population behavior around such stable states. For 
bi- and multistable systems, however, cell-to-cell variability may to a large extent be explained in terms 
of different states being occupied by different cells (even though they are genetically identical). From 
such variation different cell fates may be differentially accessible and hence understanding the causes of 
multistability in signal transaction will have ramifications across many areas of modern biology, notably 
stem-cell biology and regenerative medicine. 

One canonical class of biological systems exhibiting multistability are protein kinase cascades that 
involve multiple phosphorylation of a substrate |6|. Mitogen activated protein kinase (MAPK) cas- 
cades are the most popular exponents of this type of systemQ. Depending on the mechanism of (de- 
)phosphorylation bistability in such systems can arise [8|, which would give such systems the ability to 
use different levels of phosphorylation of the final substrate, e.g., Erk. This has been an area of growing 
activity in recent years, because of the important role MAPKs play in cell-fate determination. 

The ultimate function of Erk is to initiate a host of transcriptional responses. To fulfill such a func- 
tion, activated Erk has to shuttle into the nucleus and a growing body of recent work is paying attention 
to such spatial aspects of signal transduction ll9l4TT1l. Here we show that this spatial organization lfT2l of 
signal transduction processes plays a pronounced role in increasing the "computational space" available 
to cells. Interestingly, the same effect has been observed at the onset of mitosis lfT3l . Very much like the 
physical Address space in a computer processor [14j, the biological equivalent is influenced by the 
(bio-)physical organization of such systems. And here we show how the compartmentalization increases 
the number of stable states that can become simultaneously accessible, conferring greater flexibility and 
plasticity to such systems. Importantly, spatial organization can induce multistability into systems that 
otherwise would be mono-stable, as well as (sometimes considerably) increase the number of states in 
systems where the presence of multiple-phosphorylation sites would already give rise to multistable be- 
havior. 

This work complements the findings in llT5l by providing a detailed mathematical analysis focused 
on the Erk shuttling mechanisms. Identifying whether a system exhibits multistable behavior or not, is 
however, challenging. This type of behavior may be limited to small regions of parameter space and in 
high-dimensional spaces (our systems considered here have 20 and 36 parameters (reaction rate constants), 
respectively) it is likely not detectable by simulation or random search of the parameter space; instead 
other approaches are called for. Here we base our arguments on a set of generalizable heuristics that 
conclusively assert or reject a system's capacity for multistability, and which allow us to identify and 
delineate multi- and monostable regions in parameter space. 

1 Results 

We introduce and illustrate this framework for a basic building block of many signal transduction systems 
that for spatially homogeneous systems is guaranteed to be monostable. A one-site phosphorylation cycle 
— other post-translational modifications can, of course also be considered — consists of the reversible 
modification of a substrate S into its phosphorylated form S*. Phosphorylation and dephosphorylation are 
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Figure 1: Spatial signaling schematic. (A) Cells require movement of molecular species within the plasma mem- 
brane, cytoplasm, nucleus and more cell locations to turn genes on or off and ultimately induce a response. (B) One- 
site phosphorylation/dephosphorylation in two compartments, cytoplasm and nucleus. Molecular species: unphos- 
phorylated substrate (S), kinase (E), substrate-kinase complex (X), phosphorylated substrate (S*), phosphatase 
(F), phosphatase- phosphorylated substrate complex (Y) and superscript c denotes species in the cytoplasm and all 
others are in the nuclear compartment. Kinase and substrate total abundances (E tot: S to t) are globally conserved, 
while phosphatase abundances are conserved within each compartment (F to t, Ftot)- 



enzymatically catalyzed by E (kinase) and F (phosphatase), respectively, through a standard Michaelis- 
Menten mechanism involving the formation of intermediate complexes X, Y: 

E + S^X^E + S*, (1) 

&off,E 

F + S *^Y^F + S, (2) 

&off,F 

where k* denote the reaction rate constants. Endowed with mass-action kinetics, this cycle is monostable 

nana. 

Now suppose we learn that these enzymatic reactions can both occur in both the nucleus and the 
cytoplasm (Fig. 1) and we therefore include the reactions in the cytoplasm (denoted by the c superscript) 

E c + S c ^ X c 5^ E c + s <* (3) 

\)ff,E 

F c + S c* Y c F c + S c , (4) 

^off,F 

as well as shuttling reactions between cytoplasm and nucleus [j] 

Z c ^ Z, (5) 

^out,Z 

for the species Z of the one-site phosphorylation cycle that shuttle into the nucleus at rate fe^z and out at 
fcout,z- Total kinase, phosphatase and substrate abundances are constant for each compartment or globally 
depending on which species are allowed to shuttle between compartments. Here we use variation in the 



! We remain in a regime where compartmental models are appropriate and where we do not have to model diffusive motion 
using reaction-diffusion equations. This is appropriate for all cases where transport across a membrane is rate-limiting. 
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amount of active kinase to model how external stimuli are processed and use the substrate state to capture 
the effects of stimuli. 

The species abundances are modeled by a system of ordinary differential equations, which we analyze 
to determine if the system can exhibit multiple steady states (multistationarity) or not. This analysis 
employs a suite of different mathematical techniques, including the Jacobian injectivity criterion and 
recent developments from chemical reaction network theory |[T8l - [2T1l (see the Appendix for details). When 
multistationarity can occur we thereby obtain corresponding values of the rate constants and the steady 
states. Further, we delimit regions of the parameter space that contain all sets of rate constants that give 
rise to multistationarity. 



1.1 Analytic conditions for multistability 

Armed with these tools we establish that multistationarity cannot occur if only one species shuttles; if two 
species shuttle, then only the combinations {X, £*} or {£, Y} can induce multistationarity for certain 
total amounts and rates; and increasing the number of species that shuttle maintains multistationarity (see 
the Appendix for a full description of the sets of shuttling species inducing multistationarity). The fact 
that the model includes the formation of at least one of the intermediate complexes X, Y (and hence some 
form of sequestration) is critical for the creation of multistationarity. 

In particular, multistationarity occurs if the species E, 5, X, £* are allowed to shuttle (Fig. IB). We 
choose to study this system in detail because it corresponds to a spatial model of a simplified one-site 
MAPK model system EE) that is strictly monostable. 

If the shuttling rates fulfill 




Figure 2: Bistability in one-cycle localization model. Rate constants and total amounts are given in the Appendix. 
(A) Steady state curve of phosphorylated nuclear substrate (S*) shows bistability and hysteresis as a function of 
stimulus (total kinase E tot ). Stable states, solid lines; unstable state, dashed lines. Activation (phosphorylation) and 
deactivation (dephosphorylation) switches discontinuously from one stable branch to another at different stimulus 
thresholds, corresponding to the boundary values of the bistable regime (blue region). (B) Steady state curves of S* 
as a function of stimulus (E tot ) at varying amounts of total substrate (Stot)- (C) Steady state diagram identifying the 
regions of parameter space supporting monostability (yellow, orange) or bistability (blue) as a function of the total 
amounts of kinase (E tot ) and substrate (Stot)- (D) Theoretical western blot of bistable (two bands) or monostable 
(one band) behavior depending on dose of stimulus. 
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Figure 3: Effects of shuttling in the one-site localization model. Rate constants and total amounts as in Fig. 2. 
(A) Three steady state curve behaviors of phosphorylated nuclear substrate (S*) as shuttling are varied: high to 
low reversible bistability (blue, fc ut,E)» l° w to high reversible bistability (green, fci n ,x)> high to low irreversible 
bistability (black, fci n ,s*)- (B) Steady state curves of phosphorylated nuclear substrate (S*) as a function of the 
rate constant of S* shuttling into the nucleus (fci n ,s*) at varying amounts of stimulus (E tot ). (C) Steady state 
diagram identifying the regions of parameter space supporting monostability (yellow, orange) or bistability (blue) 
as a function of fc ou t,s an d &out,s* • 



then multistationarity cannot occur even for the spatial model, whatever total amounts and reaction con- 
stants within each cycle (see the Appendix). The shuttling rates are paired for species E, X and species 
5, and therefore a necessary condition for multistationarity is that either X moves in or out of the 
nucleus slower than the kinase, E\ or alternatively that S shuttles more slowly than its active form, 5*. 

We next home in on a set of biologically plausible rate constants and total abundances for which the 
system has three steady states, two of which are stable. We find that at low and high stimulus doses (i.e., 
total kinase level E tot = E + E c + X + X c ) there is one stable steady state of the phosphorylated nuclear 
substrate, £*, whereas for intermediate stimulus doses two stable steady states coexist (Figs. 2A-C). As 
the stimulus level changes (due to kinase production or degradation), the state of may switch either to a 
highly or lowly phosphorylated (activity) steady state based on the system's memory: switching between 
states can occur in the form of a hysteresis loop (see black and red arrows in Fig. 2A). 

Shuttling plays a pronounced role in modulating the number of discrete states of the nuclear sub- 
strate concentration that can be realized: an increase of the shuttling rate constant causes the sys- 
tem to change either from a high to a low stable state (for fc ut,E 5 ^out,x 5 ^in,s 5 ^in,s*) or vice versa (for 
fcin,E 5 Mn,x 5 ^out,s 5 ^out,s*X with an unstable steady state in between (Fig. 3 A) and saddle node bifurca- 
tions. Depending on which rate constants are altered the resulting switches can be reversible (Fig. 3A, 
fcout,E 5 &in,x) or irreversible (Fig. 3 A, fci n ,s* )• The response curves of the rate constants fci n ,E 3 ^in,s* > &out,x> 
^out,s can be tuned by altering e.g., the total kinase levels, in order to alter the size of multistabile regimes 
(Fig. 3B) or to turn reversible into irreversible switches. For example, if the shuttling rate constant fci n ,s* 
is small, then the nuclear can exist in either a high or low state (bistable region); however, following 
an increase of the shuttling rate constant across a threshold, the system switches to a monostable low 
state and cannot switch back to the high state nor re-enter the bistable regime. Similarly the rate con- 
stants fc ou t,s 5 ^in,E provide irreversible switches favoring the high state. As the shuttling rate constants 
are controlled by a variety of other processes this endows spatially structured systems with a high level of 
flexibility and increased information-processing ability compared to spatially homogeneous systems: for 
example, any violation of the inequalities, Eqns. ([6]), may result in induction of multiple stable states. 

1.2 MAPK two-site phosphorylation 

Having established that species shuttling introduces multistability in a system that otherwise is monos- 
table, we next explore the influence of shuttling in a system that can already exhibit multistability. Specif- 
ically, we consider nuclear localization in a two-site phosphorylation cycle (Fig. 4A), such as the layers of 
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Figure 4: Two site phosphorylation localization model. Rate constants and total amounts provided in the Appendix. 
(A) Schematic of two-site phosphorylation/dephosphorylation in two compartments, cytoplasm and nucleus. (B) 
Steady state curve of phosphorylated nuclear substrate (S2) shows multi stability and hysteresis as a function of 
stimulus (total kinase E tot ). Stable states, solid lines; unstable states, dashed lines. Activation (phosphorylation) 
and deactivation (dephosphorylation) switches discontinuously from one stable branch to another at different stim- 
ulus thresholds, corresponding to the boundary values of the multistable regime (blue interval). Extreme values of 
stimulus restrict the system to monostability (log scale). Closer inspection reveals up to four stable states simul- 
taneously (blue figure, zoomed, linear scale). (C, D) Steady state curve of phosphorylated nuclear substrate (S2) 
as a function of the rate constant of Si shuttling out of the nucleus (fc ou t,Si) (gray) or as a function of the total 
cytoplasmic phosphatase (F£ ot ) (coral). 



canonical MAPK cascades, and its impact on the cell's ability to establish distinct stable states. Multista- 
bility, known to exist in these systems [23] with up to three stable states ll24l . has been discussed without 
reference to any of the spatial models of MAPK signaling E El [22l [25l [26). Here, in order to differ- 
entiate between biochemical dependent and shuttling-dependent multistability, we consider biochemical 
parameter sets that preclude multistability in the absence of localization (see Appendix). Again spatial 
structure and shuttling between compartments can induce bistability; e.g., fixing the shuttling species to 
be E, Xl, X2, So, Si, S2, bistability is introduced in the system (see Appendix). Provided that the rate 
constants in at least one of the cycles are in the range of multistability of the two-site phosphorylation 
cycle, we observe that shuttling creates up to 4 stable states and 3 unstable states (Fig. 4B). Steady state 
analysis on a choice of biologically plausible shuttling rates indicates that with shuttling the two-site phos- 
phorylation cycle can undergo hysteresis effects with a large region of multistability (32 < E tot < 445), 
most of which is bistable. The doubly-phosphorylated substrate in the nucleus (S2) appears in a low/high 
monostable state for extremely low/high levels of kinase E to t (red and black arrows in Fig. 4B). 

The extended region of multistability with four stable states may provide an explanation for the versa- 
tility of MAPK signaling systems and their widespread use as relays in many signal transduction networks; 
it furthermore confers simultaneously both increased robustness and flexibility to the signal processing ca- 
pabilities of the system compared to spatially homogeneous alternatives. As expected, the steady states of 
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the system can be regulated through reversible switches governed by shuttling parameters and other total 
amounts (Figs. 4C, D). But in some cases, stable states may be so close together to be virtually indis- 
tinguishable under some physiological conditions; the switching between 3 and 4 stable states in a small 
region (see zoomed box Fig. 4B) may be an example of this. This, too, can be modulated, however, by 
regulation of the shuttling process, or by adjusting total substrate abundances. 



2 Discussion 



Fidelity of information processing and the computational capacity, i.e., the ability to map environmental 
states onto discernible internal states (in particular of proteins active in the nucleus), are thus profoundly 
affected by the spatial structure of eukaryotic cells. In a biological context, a high and stable state of 
nuclear substrate (5*) can marshal robust responses to environmental cues. To add further flexibility to 
such information processing the shuttling speed of a substrate may further depend on the nature of the 
stimulus; biological examples abound, and, for example, stimulated NIH 3T3 cells shuttle MAPK into the 
nucleus three times faster than starved cells E71l . By simultaneously controlling the stimulus dose and 
shuttling speed, a nuanced transition to reversible bistability permits hysteresis, and thus introduces the 
possibility for switching between low and high stable states (Fig. 3B). 

Trafficking is therefore more intimately related to cellular computation than is typically acknowledged 
and differences between cell lines in the shuttling rate constants of different substrates — e.g., in NIH 3T3 
mouse cells, phospho-MAPK can accumulate in the nucleus ll27ll : whereas nuclear accumulation does not 
occur in PC 12 cells [ 28 1 — may be hard- wired. Alternatively, the rich set of mechanisms affecting both 
phosphorylation (exemplary perhaps also for other post-translational modifications) and trafficking give 
cells the flexibility to change their dynamical regime, e.g., from monostable to multistable, on the fly 
and in response to further environmental and physiological cues. The spatial/compartmental organization 
of cells thus drives crucial aspects of their information processing capacity; notably the number and 
robustness of states that can be stably represented is higher (or at least as high) for spatially structured 
systems compared to homogeneous systems [lj. This provides a further rationale for the evolution of 
cellular compartments |[29l but also begs the question as to how bacteria and archaea can increase their 
computational capacities. Here, we believe, micro-environments generated by molecular crowding QUI 
confer some of the same advantages. Crowding, e.g., around membrane associated histidine kinases of 
two-component signaling systems, can isolate signal transduction components spatially from one another, 
which would have similar (although likely weaker) impact on the computational address space as cellular 
compartment have in eukaryotes. 
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A Appendix: Supporting information 

In the Supporting Information we illustrate the claims made in the main text in more detail. We further 
expand on the heuristic methods that we have developed to determine if a system of biochemical reactions 
has the capacity for multiple steady states and to find conditions on the rate constants that ensure that 



multiple steady states cannot occur. In Section A.l| we give an overview of the methods employed. In 



Section |A^2| we study a one-site phosphorylation cycle, which is monostationary, and show that shuttling 
species can introduce multistationarity. In Section |A3] we study the extended two-site phosphorylation 
cycle. Without compartmentalization the two-site modification cycle exhibits multistationarity for some 
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choices of rate constants but not all. We show that compartmentalization can introduce multistationarity 
even if the rate contents do not allow multistationarity in a two-site system without compartmentalization. 

A.l Methods for the determination and preclusion of multistationarity 

We derive conditions on the rate constants that ensure multistationarity cannot occur. Let 

/ = (/!,...,/„): R n ^R n 

be a differentiable function. The Jacobian, J x (f), of / at a point x = (xi, . . . , x n ) in W 1 is the n x n 
matrix with entry (z, j) being If / is a polynomial function then all entries of the matrix J x (f) are 
polynomials in x and consequently, the determinant of J x (f) is a polynomial in x too. The Jacobian 
injectivity criterion states that if all components fi have total degree at most two, and the determinant of 
the Jacobian does not vanish in a convex domain ft C R n , then / is an injective function in JT8ll . 

The positive steady states of a system of biochemical reactions are given as the positive solutions 
to a system of polynomial equations f K ,i(x) = 0, i = 1, . . . , n, where the coefficients of polynomials 
f K) i depend on the rate constants k = {k r } and k r is the rate constant of reaction r. If the function 
/« = > • • • 3 /«,n) is injective over the positive real numbers R™ then there can at most be one positive 
solution to f K (x) = (0, . . . , 0), and consequently, multistationarity cannot occur. We will use the Jacobian 
injectivity criterion with = R™ to determine if the function f K is injective. 

In our case, the polynomials f Kj i are either quadratic or linear, and hence the Jacobian injectivity 
criterion can be applied. The determinant of the Jacobian of f K is a polynomial in x with coefficients 
depending on the rate constants. Each term in the polynomial is of the form a^x™ 1 • • • x™ n , where a{n) 
is a coefficient and mj is 0, 1 or 2. If the coefficients of the determinant of the Jacobian are all positive or 
all negative for a specific choice of rate constants, then the determinant cannot vanish when x is positive. 
Hence, it follows from the Jacobian injectivity criterion that multiple positive steady states cannot occur. 
Importantly the coefficients a(k) themselves are polynomials in the rate constants n = {k r }. Thus, we 
can study how the signs of the coefficients vary when the rate constants vary. Our focus is on analyzing 
the coefficients in order to understand what combinations of rate constants make them all positive or all 
negative. 

The role of the Jacobian injectivity criterion is to preclude multistationarity. However, failure of the 
Jacobian injectivity criterion is not sufficient to conclude that multistationarity occurs. To investigate 
whether multistationarity occurs when the Jacobian injectivity criterion fails, we make use of the algo- 
rithms implemented in the chemical reaction network theory (CRNT) toolbox EOl . For some systems 
modeled with mass-action kinetics (as is the case here), the toolbox can conclusively determine if mul- 
tistationarity can occur or not. If the system admits multiple steady states, the toolbox outputs a unique 
set of rate constants for which there exists a pair of positive steady states (fulfilling the conservation laws 
with the same total amounts). However, the rate constants that the toolbox outputs cannot be constrained 
or controlled in any way. That is to say, we cannot ask the toolbox to restrict the search to certain re- 
gions that are considered biologically realistic. This limits substantially the use of the algorithms of the 
toolbox. Nevertheless, the output rate constants serve here as a starting point for further investigation and 
manipulation. 

In this work, a steady state is considered stable if it is asymptotically stable, that is, if the real part 
of the eigenvalues of the Jacobian of the system at the steady state are all negative. Asymptotic stability 
ensures that if the initial state of the system is sufficiently close to the steady state then it will eventually 
be attracted towards the steady state. 
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Figure 5: Shuttling of a one-site phosphorylation cycle between the nucleus and the cytoplasm. 
A.2 Shuttling in a one-site phosphorylation cycle 

Reactions and rate constants. We consider a one- site phosphorylation cycle with species 5, (the 
unphosphorylated and phosphorylated substrates), E (kinase), F (phosphatase), and X, Y (intermediate 
complexes). Phosphorylation and dephosphorylation are assumed to follow a Michaelis-Menten mecha- 
nism (see below and main text). This motif cannot admit multiple steady states and is monostable lTT6l . 



To study the effect of compartmentalization we assume that the species 5, 5*, E, X can shuttle be- 
tween the cytoplasm and the nucleus (see Figure [5}. We let Z c denote the species Z in the cytoplasm. 
Then, we have the following reactions: 



• Reactions in the nucleus: 



E + S X E + S* 

k 2 



• Reactions in the cytoplasm: 

k 7 



E c + S c 



^X c 



^E c + S c 



^4 Ate 

F + S* ^ ^ Y F + S 

k 5 



k±o fcio 
pc _|_ gc* _ u > yc ^\ pc _^ gc* 



• Shuttling reactions: 

E^^E C 

ki7 



k±4 

x- — ^x c 



k±5 



kie 

s* < s c * 

k20 



To ease the notation, we have changed the notation of the reaction constants k r in the main text 
and simply labeled them with consecutive numbers fci, . . . , k 2 o- The correspondence between the two 
notations is shown below: 



Here 


ki 


k 2 


ks 


&4 




fc 6 


A: 7 


k 8 


fc 9 


fcio 


Main text 


^on,E 




&cat,E 


^on,F 


&off,F 


&cat,F 


^on,E 


%ff,E 


^cat,E 


^on,F 


Here 


fell 


&12 


kl3 


k u 


&15 


^16 


kl7 


^18 


&19 




Main text 


%ff,F 


^cat,F 


k ou t,E 


k ut,x 


&out,S 


&out,S* 




&in,X 


&in,S 


&in,S* 



Mass-action system of ordinary differential equations. We order the set of species in the following 
way: 

(xi,x 2 ,x 3 ,X4,X5,x 6 ) = (E,X,S,S*,F,Y), (x 7 ,x$, x 9 , x 10 , x u ,x 12 ) = (E c ,X c ,S C iS c *,F c ,Y c ). 



By assuming the law of mass-action, the dynamics of this reaction network is modeled by the following 
system of ordinary differential equations (reference to time t is omitted, X{ = xi(t))\ 

x\ = —kisxi + k 2 x 2 + kzx 2 - kix\xz + /C17X7, 
X2 = -k 2 x 2 - ksx 2 - kux 2 + hxixs + hsxs, 
x 3 = k 2 x 2 - ki$xs - k±XiX3 + k 6 x 6 + fcigXg, 
x 4 = ksx 2 - h 6 X4 - /C4X4X5 + k$x 6 + foo^io? 

A 6 = /C4X4X5 - k$x 6 - k 6 x 6l (7) 

%7 = &13#1 — /C17X7 + /C8X8 + — /C7X7X9, 

X8 = fcl4^2 - ^8^8 - kgXs ~ ^18^8 + /C7X7X9, 

Xg = /C15X3 + fcsXs - ^19^9 - kjXjXg + fcl2#12 5 
*10 = ^16^4 + &9^8 - ^20^10 - fcioXioXn + k\\X\2, 
X\i = — fclo^lO^ll + ^11^12 + fcl2#12 5 
*12 = fclO^lO^ll - ^11^12 - ^12^12- 

This dynamical system has four conservation laws, accounting for the fact that the amounts of enzymes 
and substrate are conserved: 

= ±1 + ±2 + x 7 + x$, 

= x 5 + x 6 , (8) 
= ±2 + Xs + X4 + X 6 + *8 + x 9 + iio + £12, 
= i n + x\2. 

These conservation laws can be verified but adding the corresponding equations in ([7]). Since the model 
does not incorporate shuttling of the phosphatase, the amount of phosphatase is conserved separately in 
each compartment. 

Let Stot denote the total amount of substrate, and E tot , F to t, F£ ot denote the total amounts of kinase 
and phosphatase in the system, respectively. The differential equations in ([8} lead to the following equa- 
tions that are fulfilled at any time: 

Efot = Xi + X2 + x 7 + x 8 

Ftot = x$ + xq (9) 

Stot = X2 + X 3 + X4 + X 6 + X 8 + X 9 + Xio + Xu 
F tot = X U +X 12 . 

The steady states of the system are obtained by setting all derivatives ±i to zero. The system has 
the capacity for multiple steady states if there exist rate constants fci, . . . , /C20 and positive total amounts 
Stot, Etou Ftot, Ftot suc h that the equations ±i = together with ^ have more than one positive solution. 
Therefore, for fixed reaction rates and total amounts, determination of multistationarity implies solving 
a system of polynomial equations in 12 indeterminates (variables). The equations corresponding to the 
conservation laws are linear, while those corresponding to setting the derivatives to zero are quadratic 
(that is, they have terms of total degree 1 and 2). 
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Rate constants and total amounts for multistationarity (for Figs. 2, 3 in the main text). The CRNT 
toolbox provides a unique set of rate constants for which the system admits multiple positive steady states 

fei = 11.679195 k 2 = 144.94137 k 3 = 91.527059 fc 4 = 207.26904 k 5 = 22.115015, 

k 6 = 309.97808, k 7 = 49.545796, fc 8 = 8.8750284, k 9 = 262.90818, k 10 = 356.03934, 

ku = 1.8978202, k 12 = 44.457164, fc 13 = 1.0903408, k 14 = 305.42214, fc 15 = 47.547732, 

fei 6 = 41.866754, k 17 = 86.473107, fci 8 = 215.67801, k 19 = 1, k 20 = 165.98446. 

For this set of rate constants, two steady states are provided with total amounts: 

Etot = 20.7066814, S tot = 35.21053215, F tot = 3.84921092, F t c ot = 11.0903086. 

We aim to exemplify multistationarity with rate constants that are more biologically reasonable and 
of the order of experimentally determined values [9} El- To this end, we have manually investigated the 
effect of changing a specific rate or a total amount with respect to the emergence of multistationarity. We 
guide the proposed changes by the structure of the steady-state equations. This procedure has allowed 
us to tune the rate constants and total amounts to reasonable values without loosing multistationarity. 
Specifically, we settled for the rates (used to create Figures 2 and 3 in the main text): 

fei = 0.049 k 2 = 0.009 k 3 = 0.262 fc 4 = 0.356 k 5 = 0.002 k 6 = 0.044 k 7 = 0.011 
kg = 0.144 k 9 = 0.091 feio = 0.207 k n = 0.022 k 12 = 0.309 k 13 = 0.16 k 14 = 0.14 
fei 5 = 0.001 kie = 0.166 k 17 = 0.0006 k 18 = 0.33 k 19 = 0.047 k 20 = 0.041, 

and the total amounts {E to t, Stot, Ftot, Ftot} = {^2, 35, 11, 3}. With these parameters, there are three 
steady states, of which two are stable. Specifically, the steady states are approximately: 

551 = (0.676, 2.357, 16.33, 1.261, 1.023, 9.977, 17.671, 1.296, 2.068, 0.751, 2.041, 0.959) 

55 2 = (0.183, 0.965, 27.220, 0.126, 5.623, 5.377, 20.391, 0.461, 0.559, 0.107, 2.812, 0.188) 

55 3 = (1.062, 2.87, 11.847, 2.145, 0.625, 10.375, 16.371, 1.701, 3.109, 1.503, 1.546, 1.454). 

The steady state SSi is unstable and has only one eigenvalue with positive real part. 

Conditions for monostationarity (Eqns (6) in the main text). Not all choices of rate constants and total 
amounts have the capacity for multistationarity. We show here that there is a set of necessary conditions 
for the existence of multistationarity that depends exclusively on the shuttling rates. To see this, we 
apply the Jacobian injectivity criterion to a function that in part consists of the right hand sides of the 
conservation laws ([7]). 

Specifically, we consider the polynomial function f K : R 12 — » M 12 given by the right-hand side of the 
four conservation equations ([9]) and the equations in Q for all ±i except for x\, £2, £5 and in. The latter 
equations are redundant and can be obtained from the conserved equations in ([8]). The 12 components of 
the function f K = (f K l , . . . , f K l2 ) are 





= Xi + x 2 + x 7 + 


x lu 






/«,2 


= X 2 + X 3 + X4 + 


Xq + #8 


+ Xg + 


^10 + a?i2, 


fhZ,3 


= X 5 +£ 6 , 




















= k 2 x 2 - k 15 x 3 - 








Ik,6 


= ksx 2 — kiQX^ - 


k^x^x^ 


+ &5^6 " 




Ik,7 


= — k^XQ 


- k 6 x 6 , 






fhc,8 


= k 13 X! - k 17 x 7 - 


f ^8^8 ■+ 


- ^9^8 - 


k 7 x 7 x 9 , 


fn,9 


= k 14 x 2 - k$x 8 - 


/c 9 x 8 - 


k 18 x 8 + 


k 7 x 7 x 9 , 


A, 10 


k 15 x 3 + k 8 x 8 - 


/C19X9 - 


- k 7 X 7 Xg 


+ k\2X\2, 


/«, 11 


= k 16 X 4 + k 9 X 8 - 


^20^10 ■ 


- fcio^lO^ll + ^11^12 


/«, 12 


= hox^xn - knx 12 - k 12 x 12 . 
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If this function is injective over the real positive numbers WL, then multiple positive steady states with the 



same total amounts cannot occur. As described in Section \AA\ we use the Jacobian injectivity criterion 
to investigate conditions on the rate constants for which the function is injective. Since f K is quadratic, 
the criterion applies. The determinant of the Jacobian of f K can be computed using any software that 
enables algebraic (symbolic) computations, like Mathematica or Maple. We compute the determinant and 
extract the coefficients. These coefficients are polynomials in the rate constants and most of them contain 
only positive summands. Therefore, we search for the coefficients that have negative summands. After 
appropriate factorization and simplification, we conclude that the coefficients are all positive if and only 
if the following expressions are positive: 

d =k 9 k 14 + k 9 k 17 + k 3 (k 18 - k 17 ) = k 9 k 14 + (kg - k 3 )k 17 + k 3 k 18 , 

C 2 =k 3 k 12 (k 15 - k 16 )(k 18 - k 17 ) + k 3 k 15 k 16 (k 18 - k 17 ) + k 12 k 15 k 16 k 18 + k 9 k 14 k 15 k 16 

+ + k 12 k 1A ki Q k 17 + k 9 k 16 k 16 k 17 + k 12 k 1Q ki 7 ki 8 , 

C 3 =k 3 k 12 k 15 (k 18 - kir) + k 6 kgk 14 k 15 + k 6 k 12 k 14 k 15 + k 6 k 12 k 14 k 17 + k 6 kgk 15 k 17 

+ k 6 k 12 k 15 k 18 + k 6 k 12 k 17 k 18 , 
C 4 =k 3 k 15 (k 18 - k 17 )k 20 + k 6 kgk 14 k 15 + k 6 kgk 15 k 17 + k 6 kgk 14 k 20 + k 6 k 14 k 15 k 20 

+ kgki4ki 5 k 20 + k 6 kgki 7 k 20 + k 6 ki 4 ki 7 k 20 + kgk lb ki 7 k 2Q + k 6 ki 5 ki 8 k 20 + k 6 ki 7 ki 8 k 20 , 

C 5 = &3&13 + ^9(^14 - &13) + &3&18 = (&3 - ^9)^13 + &9&14 + &3&18, 

C 6 =k 6 kg(k 14 - k 13 )k 19 + k 6 k 12 k 13 k 14 + k 6 k 12 k 13 k 18 + k 3 k 12 k 13 k 19 + k 6 k 12 k 14 k 19 

+ k 3 ki 2 k 18 k 19 + k Q k 12 ki 8 ki 9l 
C 7 =kg(k 14 - k 13 )k 16 k 19 + k 3 k 12 k 13 k 16 + k 12 k 13 k 14 k 16 + k 3 k 12 k 16 k 18 + k 12 k 13 k 16 k 18 

+ k 3 k 12 k 13 k 19 + k 3 k 13 k 16 k 19 + k 12 k 14 k 16 k 19 + k 3 k 12 k 18 k 19 + k 3 k 16 k 18 k 19 + k 12 k 16 k 18 k 19 , 
C 8 =k 6 k 9 (k 14 - k 13 )(k 19 - k 20 ) + fc 9 (&i4 - k 13 )k 19 k 20 + k 6 k 13 k 14 k 20 + k 6 k 13 k 18 k 20 

+ k 3 k 13 ki 9 k 20 + k 6 ki 4 ki 9 k 20 + k 3 k 18 k 19 k 2Q + k 6 ki 8 ki 9 k 20 . 

Observe that these expressions cw/)/ involve the 8 shuttling rates and k 3 ,k§,kg, k\ 2 . Instances for which 
the coefficients Ci , . . . , Cs are negative exist. If 

&20 < ^19, &18 > &17, ^16 < ^14 > (10) 

then (7^ > for all i and hence multistationarity cannot occur for any choice of total amounts (these 
inequalities correspond to Eqns. (6) in the main text). However, there is no guarantee that when these con- 
ditions fail, the system admits multiple steady states for some total amounts. Figure [5] above is reproduced 
as Figure [6] with the shuttling rate constants indicated. 

We assume now that the dissociation constants are the same in the two compartments (the nucleus and 
the cytoplasm), that is, we assume that £3 = fcg and fcg = fci2. In this case Q > for all i 7^ 2, 8 and all 
shuttling rate constants. Therefore, two conditions suffice to guarantee monostationarity, namely: 

C 2 =fc 9 /ci 2 (£;i5 - fci6)(fci8 - fcrr) + k 9 k u ki 5 ki 6 + h 2 kuki b ki 6 

+ ki 2 k u ki 6 ki 7 + k 9 ki 5 ki 6 ki 8 + ^12^15^16^18 + ^12^16^17^18 > 0, 

C 8 =fc 9 /ci2(fcl4 - fcl3)(fcl9 _ fc 20) + ^12^13^14^20 + ^12^13^18^20 

+ k 9 kuki 9 k 2 o + h 2 kuki 9 k 2 o + k 9 ki 8 ki 9 k 2 o + ki 2 ki 8 ki 9 k 20 > 0. 
The first corresponds to C2 and the second to Cs. By inspection of these two expressions, we conclude 



that multistationarity cannot occur in any of the followin 


g cases: 






(i) 


&20 < fciQ, 


fcl8 > fcl7, 


^16 < fci5? 


fcl4 


> fe 


(ii) 


&20 > ^19, 


fci 8 > fcir, 


^16 < fc 15, 


fci4 


< fe 


(iii) 


&20 < &19, 


fcis < fcir, 


^16 > ^15, 


fci4 


> fe 


(iv) 


&20 > &19, 


fcis < fcl7, 


^16 > fcl5 3 


fcl4 


< fe 
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Figure 6: Shuttling rates for the one-site phosphorylation cycle 

Note that these only involve the rate constants for the shuttling reactions. If the dissociation rate constants 
are not exactly the same in the cytoplasm and in the nucleus, but very similar, then the conditions above 
are still sufficient. 

We see that the rate constants go in pairs: the shuttling rate constants of S relate to those of S* 9 and 
the shuttling rate constants of E to those of X. In particular, the following conditions are necessary for 
multistationarity: 

(1) If X shuttles into the nucleus slower than E then S shuttles into the cytoplasm slower than £* and 
vice versa. 

(2) If X shuttles into the cytoplasm slower than E then S shuttles into the nucleus slower than and 
vice versa. 

Sets of rate constants for which I\ < can for instance be obtained by letting the product kgkiQ be large 
and the remaining products be small such that kuki$ + ky^kyj — kukis + fcis^is < is satisfied. 



# 


No multistationarity 


Multistationarity 


1 


All 


None 


2 


{S,S*}{E,Y} {F,X} {S*,E} 
{S,F} {S,X} {S*,Y} {E,X} 
{F,Y} {S,E} {S*,F} 


{E,F} {X,Y} {S*,X} {S,Y} 


3 


{X,E,F} {Y,E,F} {X,Y,E} 
{X,Y,F} {S,F,X}{S*,E,Y} 
{S,E,F} {S*,F,E} 


{S,E,X} {S*,F,Y} {S,E,Y} {S*,E,X} 
{S*,F,X} {S, E, S*} {S*,F,S} {S, S*,X} 
{S,X,Y} {S*,Y,X} {S,F,Y} {S,S*,Y} 


4 


{Y,X,E,F} 


{S, S*,X, F} {S, S*, Y, E} {S, E, X, Y} {S*, F, X, Y} 
{S, F, X, Y} {S*,E, X, Y} {S, S*,X, Y} {S, S*,E, F} 
{S, E, F, X} {S*,E, F, Y} {S, E, F, Y} 
{S*,E,F,X} {S,S*,X,E} {S,S*,Y,F} 


5,6 


None 


All 



Table 1: One-site phosphorylation system. For all possible sets of shuttling species it is indicated if the 
system has the capacity for multiple steady states or not. 

Sets of shuttling species and multistationarity. We have shown that if the species E,X,S, 5* shuttle 
between compartments, multistationarity is created. We next investigate what the sets of shuttling species 
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Figure 7: The response plotted against the rate constant kiQ, at baseline total amounts. 

that provide multistationarity are. The results are summarized in Table [T] We use a systematic way to 
classify each motif: First, we check if the system fulfills the Jacobian injectivity criterion for all rate 
constants. If the coefficients of the polynomial in x given by the determinant of the Jacobian (as above) 
are all positive, then the system cannot exhibit multistationarity for any set of total amounts (see also 
||2T10 . If the criterion fails then we use the CRNT toolbox. 

We have obtained that if only one species shuttles then multistationarity cannot occur. That is, at 
least two species, e.g. {£*, X} or {£, Y}, are required to create multistationarity in the one-site phos- 
phorylation cycle for certain total amounts and rate constants. The addition of shuttling species maintains 
multistationarity. 

Effects of varying the shuttling rates. We analyze the steady- state response of in the nucleus as the 
shuttling rate constants and total amounts change in the system. 

The following table summarizes the type of saddle-node bifurcation curves obtained as shuttling rate 
constants of molecular species are varied. 



Rate constant 


Rate-response curve 


&13, &14, &19,&20 
^15,^17,^18 

he 


For large rate constant, only a low stable steady state is obtained 
For a small rate constant, only a high stable steady state is obtained 
Similar to the previous case, but the high branch decreases (Fig. [7]). 



By varying a total amount and shuttling rate constant simultaneously, the system may undergo irre- 
versible switches. This occurs with respect to shuttling rate constants /C14, £15, ^17, and A^o- Specifically, 
for shuttling parameters, ki$ and /C17, the irreversible switch is obtained by either increasing F tot or 
decreasing E tot ,S tot or Ff f. As the value of the shuttling rate constant increases, the response curve 
switches from a low to a high steady state, favoring accumulation in the nucleus. Conversely, increasing 
the value of shuttling parameters ku and /C20 induces an irreversible switch from the high to low steady- 
state by either decreasing F tot or increasing E to t, Stot> As highlighted in the main text (see Figure 3), the 
/C20 bifurcation is irreversible at baseline parameter values. 

A.3 Shuttling in a two-site phosphorylation cycle 

In eukaryotes, most protein phosphorylation events take place in more than one site. It is well known that 
multisite phosphorylation can cause multistationarity by itself (H[23l. However, multistationarity does 
not occur for all choices of rate constants. 

We next investigate the effect of adding species compartmentalization in a two-site (sequential) phos- 
phorylation system. We first determine rate constants for which the two-site system cannot exhibit mul- 
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tistationarity. Then, we add species shuttling and determine shuttling rate constants that induce multista- 
tionarity. 

Conditions for monostationarity in a two-site phosphorylation cycle. We consider a two-site phos- 
phorylation cycle in which modifications take place sequentially. The reactions describing the system 
are: 



S + E^±X 1 ^S 1 +E 

k 2 



^4 ha 

Sr + E X 2 ^+ S 2 + E 

k 5 



Si + F - — t Yi So + F S 2 + F ~ — t Y 2 Si + F 

kg kn 

The set of species is ordered such that 

(xi,...,x 9 ) = (E,X U S ,S U F,Y U S 2 ,X 2 ,Y 2 ). 

Assuming mass-action kinetics, then the differential equations describing the dynamics of the species 
concentrations are: 

x i = k 2 x 2 + /C3X2 - kixixs - /C4X1X4 + k$xg + k^xs, 
x 2 = -/c 2 x 2 - + kixixs, 
X3 = k 2 x 2 — kixixs + kgXQ , 

X4 = /C3X2 — /C4X1X4 — /C7X4X5 + fcsX6 + + ki 2 Xg, 

X5 = — /C7X4X5 + + ^9 X 6 — kioX$X7 + fcuXg + k\ 2 X$, 

Xq = /C7X4X5 — — kgXQ^ 

x 7 = -ki x 5 x 7 + k 6 x 8 + fcuXg, 
x's = /C4X1X4 — /C5X8 — kfiXg, 

Xg = /C10X5X7 - fen £9 - fci2#9- 

This system has the following conserved amounts: 

Efot = X\ + ^2 + ^8, ^ot = ^5 + Xq + £9, = X 2 + £3 + X 4 + X 6 + X 7 + X 8 + X 9 . 

The steady-state equations are given by X{ = 0. Because of the constraints given by the conservation 
laws, the equations x\ = 0, x 2 = and x$ = are redundant and can be removed. 

We proceed as above to determine rate constants for which the system cannot have multiple steady 
states. That is, we apply the Jacobian injectivity criterion. We consider the function f K : R 9 — >► R 9 given 
by the three equations coming from the conservation laws and the 6 remaining steady-state equations: 



f^i( x 

f^( x 

f*fi( x 
f*j( x 
f*A x 



= X\ + x 2 + x$, 

= X 2 + X 3 + X 4 + Xq + X 7 + X 8 + X 9 , 

= X5 + X6 + Xg, 

= k 2 X 2 — kiXiXs + kgXQ, 

= hx 2 - /C4X1X4 - /C7X4X5 + /c 8 x 6 + hxs + fci 2 Xg, 

/C7X4X5 - k 8 x 6 - k 9 x 6 , 
= -ki x$x 7 + k 6 x 8 + fcnXg, 
— /C4X1X4 — — ^6^8, 

= /C10X5X7 - fcnXg - ki 2 X 9 . 
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If this function is injective over the positive real numbers then multiple positive steady states cannot occur 
with the same total amounts in the conservation laws. Next, we compute the determinant of the Jacobian of 
f K . As a polynomial in x, the only coefficients (which depend on the rate constants ki) of the determinant 
that are not sums of positive terms are: 

C\ — —k2k±k§k 7 k§k\§ — ksk^k^fak^kio — k\k^k§kjk§kii — k\k^kQk 7 kgki2 

+ hksk^kioku + kiksk 6 k 7 kioki2 + k\k^k^k%k\^k\2 + kifak^k^k^ 
C 2 = -kik4k 7 ki (k 6 k 9 - k 3 ki 2 ). 

If a choice of rate constants fulfills C\ , C 2 > 0, then the system cannot have multiple positive steady 
states for any total amounts. The coefficient C\ can be rewritten as 

Ci = -k^k^^kio + /C3/C10 + k\kn + k\kvi) + ^3^12^1^10(^5/^7 + k 6 k 7 + k^ks + k^kg). 

For C2 > we require 

k 3 /k 6 > kg/k 12 . 
If this inequality is fulfilled then C\ > if also 

/c4/c7(/c 2 /cio + kskio + kikn + hku) < kikio(k$k 7 + k 6 k 7 + /c 4 /c 8 + ^4^9)- 

This inequality can be rewritten as 

k 2 + k 3 fcn + fcig < fc 5 + fc 6 fcs + k 9 
k\ kio k4 k 7 

Let 



If 



fa + &6 &2 + fa k 8 + fcg fcn + fa kg 

Oil — ; ; , OL 2 — , a 3 = — . 

fa ki k 7 fcio fa ku 



OL\,OL 2 , «3 > 0, (11) 



then the two-site phosphorylation system cannot have multiple positive steady states no matter the values 
of the total amounts. That is, a sufficient condition for the preclusion of multistationarity is obtained. 

Observe that ol\ and a.2 imply an inequality between the Michaelis-Menten constants of the kinase 
and the phosphatase in each phosphorylation site. Namely, the Michaelis-Menten constant of E for the 
second site is larger than that for the first phosphorylation site, and the Michaelis-Menten constant of F 
for the first site is larger than that for the second site. 

Negation of this condition is a priori not sufficient to guarantee multistationarity. However, if 0^3 < 0, 
then there exist total amounts for which the motif exhibits multistationarity 1(311 . 

Mass-action system of ordinary differential equations for the two-site shuttling. Consider now two 
copies of a two-site phosphorylation cycle as above and let So, Si, S2, X\, X 2 , E shuttle between the 
nucleus and the cytoplasm. The reactions describing the system are: 

• Reactions in the nucleus: 

S + E Xi S ± + E S ± + E X 2 S 2 + E 

k 2 k$ 

&7 ho ^10 fcio 

S\ + F zz=zz Yi S + F S 2 + F z==^ Y 2 S! + F 

kg ku 
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Reactions in the cytoplasm: 

s% + e c - — t x{ 5f + £ c 5f + £ c i — ^ x| s% + £ c 

kl9 fo. k 2 2 kod 

S{ + F c ~ — * Yf S§ + F c ^ + F c ~ — * Y 2 C S{ + F c 



&20 &23 



Shuttling reactions: 



&25 &26 &27 

/C3i /C32 """ &33 



&28 &29 &30 



&34 ' &35 &36 



The species are ordered as: 

(x 1 ,...,x 9 ) = (E,X 1 ,So,S 1 ,F,Y 1 ,S 2 ,X 2 ,Y 2 ) 

(x w , x 18 ) = (E c , XI, S c , SI, F c , Y{, S c 2 , X c 2 , Y 2 C ). 

Assuming mass-action kinetics, the system of differential equations describing the dynamics of the 
centrations of the species is: 

±i = k 2 x 2 + ksx 2 - kixixs - /C4X1X4 + /c 5 x 8 + k 6 x 8 - k 2 ^x\ + A^i^io, 
x 2 = -k 2 x 2 - k 3 x 2 + kixixs - k 26 x 2 + h&xn, 
x 3 = k 2 x 2 - kixix 3 + k 9 x 6 - k 27 x 3 + k 33 xi 2 , 

x 4 = k 3 x 2 - /C4X1X4 - £7X4X5 + fc 8 x 6 + fc 5 x 8 + h 2 xg - /c 28 x 4 + k M xi 3 , 

X5 = -/C7X4X5 + fc 8 X 6 + /c 9 X 6 - /C10X5X7 + fcuXg + fcl2^9, 
X 6 = £7X4X5 - fc 8 X 6 - k 9 X 6l 

x 7 = -/C10X5X7 + fc 6 x 8 + knxg - k 29 x 7 + k 3 $x 16l 
x 8 = k±xiX4 - fc 5 x 8 - fc 6 x 8 - k 30 x 8 + k 36 x 17l 

Xg = /C10X5X7 - fcuXg - fci2#9, 
*10 = ^14^11 + ^15^11 - ^13^10^12 - ^16^10X13 + ^7X17 + fci 8 Xi 7 + k 2 ^X\ - /C31X10, 
Xn = -kuXn - /C15X11 + ^13^10^12 + ^26^2 ~ ^32^11, 
%12 = ^14^11 - ^13^10^12 + ^21^15 + ^27^3 ~ ^33^12, 

Xl 3 = kisXn - ki 6 XioXi 3 - /C19X13X14 + k 2 ()Xi5 + /C17X17 + /C24^18 + ^ 28 X4 - £34X13, 
X14 = -^19^13^14 + ^20^15 + ^21^15 ~ ^22^14^16 + ^23^18 + ^24^18, 
Xl5 = ^19X13X14 - fc 2 0^15 - ^21^15, 

X 16 = -k 22 X U XiQ + fci 8 Xi7 + fc 2 3^18 + ^29^7 ~ ^35^16, 
Xtf = fci 6 Xi Xi3 - ki 7 Xi 7 - fci 8 Xi7 + /C 30 X 8 - /C36X17, 
*18 = ^22^14^16 - ^23^18 ~ ^24^18- 

This system has the following conserved amounts: 

E t ot = xi + x 2 + x 8 + X10 + Xn + X17, 
Ftot = x$ + x 6 + x 9 , 

^ot = ^14 + ^15 + ^18, 

Stot = X 2 + X 3 + X 4 + X 6 + X 7 + X 8 + X 9 + Xn + X12 + X i3 + X i5 + X i6 + X i7 + X i8 . 
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Creation of multistationarity in a two-site phosphorylation cycle. We use the CRNT toolbox to obtain 
initial rate constants and total amounts for which the system admits multiple positive steady states. The 
output rate constants do not fulfill 011,012,013 > in each compartment independently. Hence, it is 
not possible to decide whether multistationarity arises due to shuttling or due to phosphorylation of two 
different sites. 

Next we investigate the effect of changing the rate constants with respect to the existence of mul- 
tistationarity for the same total amounts. We proceed by manually modifying the rates while keeping 
multistationarity and such that the sufficient conditions for the preclusion of multistationarity in a two-site 
phosphorylation cycle are satisfied in each compartment. 

We end up with the following rate constants: 

• Reaction rates in the nucleus: 

fci = 100 k 2 = 2 k 3 = 10 k 4 = 80 k 5 = 6 k 6 = 6 

k 7 = 350 k 8 = 3 k 9 = 10 fcio = 650 k u = 8 k 12 = 8. 

• Reaction rates in the cytoplasm: 

k 13 = 300 k u = 1 fci 5 = 10 k 16 = 50 k 17 = 1 k 18 = 1 

k 19 = 350 k 20 = 30 k 21 = 190 k 22 = 150 k 23 = 2 k 24 = 20. 

• Shuttling rates: 

k 25 = 10 k 26 = 30 k 27 = 70 k 28 = 30 k 29 = 1 k 30 = 10 
k sl = 450 k 32 = 20 k 33 = 20 k M = 25 k 35 = 10 k 36 = 100. 

This choice of rate constants fulfills that ai, a 2 , a 3 > in the nucleus and in the cytoplasm (where in 
the later, indices of the rate constants in a* are shifted by 12). Therefore, with these rate constants, the 
two-site phosphorylation cycles in the nucleus and in the cytoplasm cannot have multiple positive steady 
states independently of each other. 

The system with the shuttling reactions, however, does have the capacity for multistationarity. Specif- 
ically, if the total amounts are set to: 

Etot = 50, Stat = 100 F tot = 15 F t c ot = 21, 

then the system has three positive steady states: two of them are stable and one is unstable. The positive 
steady states are the following: 

551 =(1.89, 9.18, 0.62, 1.88, 0.01, 0.7, 25.11, 23.21, 14.28, 0.07, 13.38, 6.8, 0.04, 18.49, 1.15, 0.01, 2.28, 1.36) 

55 2 =(6.45, 5.93, 0.1, 0.61, 0.01, 0.17, 35.4, 25.64, 14.82, 0.13, 9.33, 2.8, 0.03, 18.32, 0.8, 0.02, 2.52, 1.89) 

SS S =(0.37, 17.59, 5.96, 2.3, 0.17, 10.66, 0.6, 5.65, 4.16, 0.04, 25.8, 24.89, 0.06, 19.22, 1.72, 0.00044, 0.56, 0.06). 

Here SSi is the unstable steady state. 

Rate constants to obtain seven steady states (for Figs. 4(B-D) in the main text). Additionally, we 
have observed that, with specific choices of rate constants and total amounts, up to 7 steady states can be 
created in this system. For instance, consider the rate constants (used in the main text to create Figures 
4(B-D)): 

(fci, . . . , k 12 ) = (101, 2, 11, 79, 6, 6, 568, 3, 12, 1502, 8, 8), 
(fc 13j . . . , k 24 ) = (210, 3, 34, 49, 1, 1, 344, 26, 187, 149, 1, 1), 
(fe 2 5, • • • , ^36) = (10, 34, 7.5, 312.5, 0.1, 0.1, 44, 23, 2, 250, 0.1, 0.1), 
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Figure 8: The steady state response curve of S2 as stimulus (E to t) varies. All other parameters are at 
baseline values. Stable state, solid line; unstable state, dashed line. 

and the total amounts (Etot, Stot, E to t, F^ ot ) = (57, 111, 15, 21). This system has 7 steady states, 4 of 
which are stable: 

SS t = (18.27, 2.41, 0.04, 0.03, 5.05, 6.28, 0.01, 4.12, 3.67, 5.2, 1.45, 0.005, 0.21, 0.01, 0.002, 46.2, 25.61, 20.99) 

55 2 = (12.66, 1.2, 0.03, 0.17, 0.47, 3.11, 0.26, 14.67, 11.42, 3.4, 0.73, 0.004, 0.3, 0.01, 0.004, 33.77, 24.34, 20.99) 

55 3 = (13.01, 0.05, 0.0009, 0.23, 0.01, 0.09, 16.34, 19.9, 14.9, 2.96, 0.06, 0.003, 0.29, 0.02, 0.01, 17.14, 21.06, 20.1) 

55 4 = (13.66, 2.04, 0.004, 0.23, 0.004, 0.039, 35.84, 20.54, 14.96, 2.63, 3.95, 0.3, 0.22, 2.48, 0.84, 0.1, 14.18, 17.76) 

55 5 = (14.56, 3.53, 0.006, 0.22, 0.004, 0.03, 40.64, 20.63, 14.96, 2.49, 6.85, 0.56, 0.14, 6.54, 1.45, 0.03, 8.95, 13.01) 

55 6 = (4, 15.03, 0.4, 0.49, 0.27, 5.07, 0.38, 12.88, 9.65, 0.12, 23.78, 35.81, 0.2, 14.89, 4.88, 0.001, 1.19, 1.23) 

55 7 = (5.94, 11.14, 0.15, 0.5, 0.02, 0.43, 6.89, 19.51, 14.55, 0.12, 18.91, 31.22, 0.17, 14.86, 4.07, 0.002, 1.38, 2.07) 

The stable steady states are SSi, SSs, SS5, SSq. 

As Etot varies, the number of states changes as shown on a log-log scale in Figure [8| corresponding 
to the bifurcation diagram (Figure 4B in the main text, semi-log scale). 

B Software tools 

Calculations to assess multistationarity were made using Mathematica (Wolfram Research, Inc., Mathe- 
matica, Version 7.0, Champaign, IL., 2010) and CRNT toolbox [20]. Bifurcation diagrams were computed 
using Oscill8 E3 and visualized with MATLAB (The MathWorks Inc., Natick, MA, R201 lb). Figures in 
the main text were created using Adobe Illustrator (Versions CS3 and CS4. San Jose, California: Adobe 
Systems, Inc., 2011). 
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